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Abstract. The localization of populations of replicating bacteria, viruses or 
autocatalytic chemicals arises in various contexts, such as ecology, evolution, medicine 
or chemistry. Several deterministic mathematical models have been used to characterize 
the conditions under which localized states can form, and how they break down due 
to convective driving forces. It has been repeatedly found that populations remain 
localized unless the bias exceeds a critical threshold value, and that close to the 
transition the population is characterized by a diverging length scale. These results, 
however, have been obtained upon ignoring number fluctuations ("genetic drift"), which 
are inevitable given the discreteness of the replicating entities. Here, we study the 
localization/delocalization of a finite population in the presence of genetic drift. The 
population is modeled by a linear chain of subpopulations, or demes, which exchange 
migrants at a constant rate. Individuals in one particular deme, called "oasis" , receive 
a growth rate benefit, and the total population is regulated to have constant size N. 
In this ecological setting, we find that any finite population delocalizes on sufficiently 
long time scales. Depending on parameters, however, populations may remain localized 
for a very long time. The typical waiting time to delocalization increases exponentially 
with both population size and distance to the critical wind speed of the deterministic 
approximation. We augment these simulation results by a mathematical analysis that 
treats the reproduction and migration of individuals as branching random walks subject 
to global constraints. For a particular constraint, different from a fixed population size 
constraint, this model yields a solvable first moment equation. We find that this solvable 
model approximates very well the fixed population size model for large populations, but 
starts to deviate as population sizes are small. Nevertheless, the qualitative behavior 
of the fixed population size model is properly reproduced. In particular, the analytical 
approach allows us to map out a phase diagram of an order parameter, characterizing 
the degree of localization, as a function of the two driving parameters, inverse population 
size and wind speed. Our results may be used to extend the analysis of delocalization 
transitions to different settings, such as the viral quasi-species scenario. 



Contents 



1 Introduction [2] 

2 Simulations [5] 

2.1 Simulation model [5] 

2.2 Simulation results IH] 

3 Analytical approach 1101 

3.1 Deterministic limit [H] 

3.2 General approach to finite population sizes [TT] 

3.3 Stochastic oasis model [TU 

3.3.1 Stochastic oasis model without convection [15] 

3.3.2 Deterministic and stochastic limits [IT] 

3.4 Stochastic oasis model with convection [TH] 

3.4.1 Approximation for the death rate [19] 

4 Discussion [21] 
Acknowledgments [22] 
A Different choices for the definition "population on the oasis" 1221 
B Relaxation spectrum in the case without convection 1231 
References 1261 



1. Introduction 



Many population biological models simplify ecological complexity by assuming that 
habitats are homogeneous. Within these models, populations have a strong tendency 
to become evenly distributed in space as time proceeds. However, many biologically 
plausible habitats are strongly heterogeneous and provide gradients or localized spots of 
growth rather than homogeneous growth conditions. For instance, hydrothermal vents are 
highly localized environments that may have played a key role in the evolution of life itself 
[23]. Today chemo-synthetic bacteria using the heat, methane and sulfur of such vents 
form the basis of deep sea food chains, ranging from bacteria up to higher organisms as fish 
or octopuses [2"Tll2"2"] (Figure la). An oasis in the desert is another example, where member 
species shape the environment for one another, thereby feeding back on the stability of the 
oasis (Figure lb). When such heterogeneities are included in the population dynamics, 
the question arises whether populations will also become heterogeneously distributed, or 
localized, despite the mixing forces of diffusion or convection. 
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Figure 1: Examples for strongly heterogeneous environments, (a) Chemo- 
synthetic bacteria near hydrothermal vents. (Photo: NOAA) (b) Localized growth in 
a desert: an oasis. (Photo: Luca Galuzzi) 

There have been several theoretical accounts to model the population dynamics at 
localized growth spots [21 [201 13] • These studies focused in particular on the effect of 
a convection term, or wind, on the population localized at a growth spot, called the 
"oasis" . It has been found found that populations are localized as long as the convection 
forces remains below a certain critical value. As the convection speed is increased, the 
population becomes more widely spread around the oasis. The spread of the population is 
characterized by a correlation length scale that diverges at the critical convection velocity, 
signaling a continuous phase transition to extended states. This derealization transition 
has also been investigated experimentally taking advantage of the fast growth rate of 
microbes [TJJJ [TB|. In these experiments, bacterial colonies were grown under hazardous 
UV light. Only a small localized (and moving) plate shielded the colonies, thus creating 
a preferred growth environment similar to the oasis considered in the theoretical models. 

Closely related localization phenomena have been studied in evolutionary biology. 
In this context, populations evolve in a heterogeneous fitness landscape, rather than a 
geological landscape. Deleterious mutations act like a convection force that tends to push 
populations away from the optimal (i.e. most fit) genotype. For small mutation rates, 
natural selection is able to maintain a population cluster localized at the fitness optimum. 
For large mutation rates, however, populations drift away from the peak. The associated 
derealization transition is called genetic meltdown in the population genetic literature 
[HJ EE], respectively error threshold in the viral quasi-species literature [5j El [101 l2T| [28] . 

The results on both the ecological derealization transition and the error threshold, 
have been obtained largely by a linear stability analysis of the localized state. This 
approach ignores number fluctuations, which are however inevitable in populations of 
discrete individuals. The fluctuations are generated by variations in the offspring number, 
which itself depend on factors such as environment, lifespan or genetic background. The 
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Figure 2: Hypothesized phase diagram for finite populations localized at an 
"oasis" . Deterministic calculations of a population localized at a favorable growth spot 
("oasis") predict a sharp derealization transition: When the convection speed exceeds 
a critical value v c , populations escape from the oasis. Our stochastic simulations of 
populations of fixed size suggest that this picture has to be modified to account for 
number fluctuations, i) Populations always delocalize on sufficiently long time scales, due 
to number fluctuations. The escape time, however, depends exponentially on population 
size. Thus, derealization is effectively unobservable in sufficiently large populations, 
ii) In the presence of convection, the escape times are strongly reduced. In an exactly 
solvable model with fluctuations but finite populations, a sharp transition line can be 
derived (red dashed line). This line also describes qualitatively the crossover from small 
to very large escape times (compared to the coalescence time) in fixed population size 
models. 
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question thus arises how the standard view of the localization transition has to be modified 
to account for number fluctuations. We hypothesize that localized populations can be 
destabilized not only by convection but also by number fluctuations. Accordingly, we 
expect a modified phase diagram that depends on convection speed and population size, 
such as sketched in the tentative phase diagram (Figure [2]). 

To test our hypothesis, we study quantitatively the localization transition in the 
presence of number fluctuations. We first simulate a one-dimensional population localized 
near a single lattice site with enhanced growth rate. This oasis is surrounded by a 
homogeneous habitat of lower growth rate. In order to fix the total population size, 
the excess population produced in each generation is removed by a (space-independent) 
death rate, which is adjusted in each time step. Our main goal is to reveal the impact 
of number fluctuations on the degree of localization, and the population distribution at 
steady state. We augment our simulation results by a mathematical analysis based on 
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Figure 3: Stochastic dynamics of the "oasis" model. A single time-step in our 
simulations consists of two sub-steps, (a) Individuals migrate between neighboring lattice 
sites, or demes, in general with a bias in one direction to account for convection. Number 
fluctuations are incorporated to account for the population turn over from one generation 
to the next. The subpopulation located at the oasis enjoys an increased reproduction rate 
a. This first sub-step effectively describes a branching random walk with a branching 
rate that depends on the location, (b) In the second sub-step, the (global) population 
is regulated to comply with the constraint of fixed population size (see main text for 
details). 



branching random walks subject to a tuned global constraint. This solvable model has 
been shown to be a good approximation to fixed population size models as far as the 
description of noisy traveling waves is concerned [5]. 

2. Simulations 

2.1. Simulation model 

In our simulations, the population of reproducing individuals is distributed along a linear 
chain of sub-populations, called demes [HI [25]. These demes are considered to have a 
distance d = 1 to their neighbors. The number of individuals in deme i is denoted by q. 
Each time step of duration e (typically e = 0.2) generations consists of the several sub- 
steps illustrated in Figure [3| First, we add to each occupancy variable c$ the deterministic 
change A^ e ^Cj given by 

A det Q/e = (cj+i + Ci_i - 2ci) + v (c i+1 - C;_i) /2 + a5 i0 c . (1) 

The first two terms on the right hand side account for biased migration: Individuals jump 
at rate l + v/2 from % to % + 1 and at rate l — v/2 from % to % — 1 . The last term represents 
deterministic growth at the oasis [i = 0) at rate a. 
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Next, we account for a stochastic change by adding A 



Q given by 



A stoch Ci /v^ = V2 (Poisson(ci) - a) . 



(2) 



Here, Poisson(cj) is a random number drawn from a Poisson distribution with mean q. 
The particular form of the noise in ^ ensures that i) occupancies do not drop below 
zero and ii) the generated variance in occupancy is given by 2eCj. The variance produced 
thus equals the expected variance in a time e for a (large) population of individuals 
that duplicate and die at the rate 1 (in our computational unit of time). We note 
that the natural alternative of directly simulating individuals that duplicate and die is 
computational much more expensive if is much larger than 1. 

After the previous sub-steps, the population size will typically be larger than the pre- 
set value N. In order to restore the population size constraint, we add to each variable 
Cj the change 



The rescaling of the occupancies by the factor of 1 — r effectively represents a homogeneous 
culling of the population that is independent of the location. Notice that, although our 
simulations account for the fluctuations induced by the discreteness of the individuals, the 
occupancies Cj are generally non-integer. The approximations involved are equivalent to 
the diffusion approach widely used in population genetics [131 S]: As long as populations 
are large (N 3> 1) and the temporal changes exhibit the same mean and variance as the 
discrete particle system, the deviations between both approaches should be small in the 
limit e — > 0. This means in particular that the results become insensitive to the precise 
noise distribution used, as long as the first two moments are correct. 

We impose periodic boundary conditions, as if demes were arranged on a circle like 
in the experiments of Ref. [19J. The number of demes is chosen large enough to ensure 
that population fits the simulation box. Time averages are taken after the simulations 
have become independent of their localized initial conditions. 

2.2. Simulation results 

Snapshots of the population density are shown in Figure |4j For the chosen parameter 
values, fluctuations in the population density are clearly visible. The time traces 
of the center-of-mass trajectory (lower half of Figure [1]) indicate that the population 
peak fluctuates wildly around the oasis, until it eventually escapes. The escape time 
(or derealization time) is much greater in the case without convection (a) than with 
convection (b), but it eventually occurs. In fact, no sharp transition between localized 
and delocalized populations is detectable in our noisy system. If the population size is 
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Figure 4: Snapshots of simulations with fixed population size. Figure (a) depicts 
snapshots for our simulation model, in which individuals carry out an unbiased branching 
random walk subject to the global constraint of fixed population size (see Figure [3] for a 
description of the time step). All lattice sites (labeled by an index i) are equivalent, except 
for deme i = where the growth rate is increased by a = 0.1. The resulting density profiles 
are typically quite narrowly centered around the oasis (c.f. green and red lines). The 
average profile (black line) is much broader than typical profiles, due to large excursions 
of the center-of-mass. These excursions can be seen in the center-of-mass trajectory at the 
bottom part of the figure. Notice that the population occasionally escapes from the oasis, 
due to a rare large fluctuation. Such derealizations are recorded when the occupancy 
at the the oasis drops below a critical value, en < 1. In the convection-less case, the 
population returns quickly after a random excursion (e.g. the first derealization event 
in (a)). With convection, (b), the population usually traverses the (periodic) simulation 
box after a derealization event until it reaches the oasis again. The escape time between 
two derealizations is indicated by a black arrow. The green and red bullets indicate 
the spatio-temporal position of the center-of-mass corresponding to the red and green 
snapshots in the upper part of the figure. The convection speed in (b) is set to 30% of the 
critical convection speed v c = ad, at which derealization is observed in the deterministic 
limit [2] . Notice that convection leads to more frequent escapes from the oasis. Simulation 
parameters were a = 0.1, iV = 3000, e = 0.2. 



large enough, the derealization is however unobservable within reasonable (simulation-) 
time. For definiteness, we say the population is on the oasis, when at least one individual 
is present there, en > 1, and off the oasis, when this condition is not met. We checked 
several alternative conditions, they yield similar or comparable results (see appendix |A|) . 
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Time averaged quantities (i.e. occupancies) are then obtained by averaging only over 
states with c > 1. 

We measured in detail the time-averaged death rate p = rj (a 2 d 2 /AD), scaled by 
its value a 2 d 2 /4D in the deterministic no-convection limit. If the overall death rate is 
large, a population can only be sustained in the vicinity of the oasis where net-growth is 
positive. Consequently, populations will be strongly localized at the oasis for p = 0(1). 
For smaller and smaller death rates, the population will be able to explore larger and 
larger parts of the habitat. When the death rate eventually changes sign (thus becoming a 
growth rate), growth becomes possible throughout the entire habitat. As a consequence, 
the population distribution becomes necessarily extended in space. Therefore, the death 
rate p may be considered as an order parameter that measures the degree of localization 
of the population at the oasis. Consistent with this interpretation, the mean field theory 
predicts that p continually decays from 1 to as the convection speed is increased from 
to the critical value v c = ad. 

The death rates p obtained from simulations are depicted in Figure |5j across a range 
of population sizes and convection velocities bounded by v c = ad for which derealization 
is observed in the deterministic limit. As expected from the mean field theory, we find that 
the death rate continually decays to zero with increasing convection v. In addition, we find 
that death rates decrease when we decrease the population size. For any finite population, 
death rates go to at convection speeds significantly smaller than expected from the mean 
field predictions. Thus, increasing the demographic noise tends to destabilize the localized 
populations. 

Another more direct measure for the degree of localization is the escape time T 
between a localization and the consecutive derealization event, as indicated in the the 
center-of-mass trajectories in Figure |4j Histograms of escape times for various population 
sizes are shown in Figure |6a| Notice that the escape time distribution is to a good 
approximation exponential, 

Pr(T) ~ exp (--) , (5) 



T f 

indicating a memory-less derealization process with a characteristic escape time r. 
Deviations from the exponential distribution are discernible only on short times (inset 
in Figure 6a). These deviations result from fluctuations in the density profile that 
temporarily satisfy the derealization condition cq < 1. After a relaxation process, the 
occupancy at the oasis is restored. Genuine derealization is therefore described by only 
the exponential part of the distribution. The characteristic escape time r can be extracted 
from the slope on a semi-logarithmic plot. This quantity r is depicted in Figure 6b as 



a function of population size for various convection velocities. Notice that r decreases 
when either the population size is decreased (because the relative strength of fluctuations 
increases), or the convection speed is increased (because convection removes individuals 
from the oasis, thereby spreading the population). Numerically, our observations can be 
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Figure 5: Mean death rates for varying population size and convection speed. 

In our simulations, the death rate r is the factor with which the population has to 
be scaled down to keep its global size constant (c.f. section 2.1). Here, we display 
the scaled death rate p = f / (a 2 d 2 /AD) normalized by its deterministic no-convection 
limit a 2 d 2 /4D. According to the deterministic limit, the death rate should decay with 
increasing convection speed according to p = 1 — (v/v c ) 2 (solid black line) up to the 
critical convection speed v c = ad [2]. For finite populations the death rates are shifted 
to lower values roughly by the constant 2/N, as predicted by the analytical approach 
explained in section |3j As a consequence, smaller convection speeds are required for 
the death rates to cross 0, which signals the derealization of the population. Our time 
average only includes states in which there is at least one individual at the oasis, cq > 1. 
For convection speeds larger or equal than the deterministic critical velocity, this is just a 
tiny minority of all states. Hence, the points for large convection speeds (close to v c ) are 
slightly biased. Solid lines are obtained by our theoretical approach explained in section 
[3l in particular by the relation (48). 



summarized by (Figure 6b) 



r ~ exp I — 1 (6) 



V v 

with 77 « f3 (1 — v/v c )~ and (5 « 560. 

In summary, our numerical simulations indicate that, due to number fluctuations, 
populations inevitably escape the "pinning" at the oasis on sufficiently long time scales. 
The typical escape times, however, increase exponentially with population size, and thus 
derealization is effectively unobservable for sufficiently large populations in the absence 
of convection. Our goal for the following is to develop an analytical framework to predict 
a critical population size in order to observe localization for a given convection velocity. 
This will be used to support the phase diagram in Figure [2] 
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Figure 6: Distribution of escape times T measured in simulations of fixed 
population sizes, (a) Distribution of escape times T for different population sizes in the 
absence of convection. In all cases, the distribution is close to exponential. Deviations 
at early times (inset) are attributed to shape fluctuations, see main text. From the 
exponential tail of the distribution, a characteristic escape time r can be extracted. This 
time strongly depends on the population size and convection velocity, as shown in figure 
(b). Solid lines in figure (b) represent the fit in (|6]). The three black arrows indicate 
the characteristic times r of the distributions shown in (a). The simulation results were 
obtained using parameter values a = 0.1, e = 0.2. 

3. Analytical approach 

While our simulations are discrete in both space and time, they should be amenable to a 
continuous description if the time step e and the reaction rate a are small enough. Under 
these assumptions, the stochastic dynamics can be approximated by a reaction-diffusion 
equation of the form 



Here, D is a diffusion coefficient, which has to be set equal to 1 deme 2 per generation to 
match our simulations. The second term on the right hand side describes convection with 
velocity v. The last "reaction" term K(x, t) summarizes all stochastic and deterministic 
changes due to the turn-over of the population. An example of the successful application 
of such react ion- diffusion equations in population biology is the Fisher-Kolmogorov wave 
equation, in which the reaction term describes logistic growth [15]. This equation 
exhibits traveling wave solutions, which have been used to describe range expansions, 
epidemic outbreaks or the spreading of beneficial mutations [12], [HI [9] . 




(7) 
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3.1. Deterministic limit 



In the limit of very large populations, we can use a deterministic reaction term to describe 
the population dynamics in the oasis model, 

K(x,t) = [a8(x) -r]c(x,t) . (8) 

The first term approximates the growth rate benefit provided at the oasis by a Dirac- 
delta-function 8(x). The pre-factor a represents the product of the increased growth rate 
a at the oasis and the linear size d of the oasis, a = ad. The constant r > represents 
the death rate, which is tuned such that (]7]) admits a stable steady state. 

The analysis of the resulting reaction-diffusion system has been carried out in Ref. [2J. 
It is found that the rescaled death rate p = ADr/a 2 = 1 in the absence of convection, 
v = 0, and continually decays to as the convection speed approaches the critical speed 
v c = a, 

p(v) = 1 - {v/v c f . (9) 

It can be seen from Figure [5] that this mean field result clearly overestimates the death 
rates measured in the stochastic simulations. 

The localized steady state profile is found to be given by 

c st (x) ~ exp [(-a|x| + vx)/(2D)} . (10) 



The density profile in equation (10) has characteristic width £(i>), given by 

= — , (ii) 

a — v 

which diverges as a power law as the convection speed approaches the critical speed. This 
diverging length scale indicates that convection driven derealization is a continuous phase 
transition in the deterministic limit. 

The predicted density profiles are depicted in Figure [7j together with stochastic 
simulation data. The agreement is for the chosen parameter sets rather good, which is 
in contrast to the discrepancy between predicted and measured death rates. 



3.2. General approach to finite population sizes 

Modeling the spatial dynamics of finite populations typically leads to stochastic non- 
linear reaction diffusion equations: A non-linearity is necessary to fix the population 
size, and stochasticity arises due to the finiteness of the population. As a consequence, 
these models are generally hard to treat analytically. Systematic approximation schemes, 
such as perturbation expansions, often fail, as the example of noisy traveling waves 
demonstrates [23] . 

In a recent study, we have shown, however, that there is one particular form of the 
non-linearity that allows treating these and similar problems analytically [§]. We now 
review our methodology, and describe how the same method can be applied to diverse 
situations, including both ecological and evolutionary problems. 
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Figure 7: Density profiles obtained for different models and approximations. 

(a) and (b) correspond to cases without and with convection, respectively. Symbols 
represent simulations with population size constraint. Green lines are predictions for the 
tuned constraint model, which we can derive analytically in the convection-less case (a). 



When including convection, only approximations are possible (b) (c.f. section 3.4). Blue 
lines represent the deterministic mean field approximations. For the parameters used 
(a = 0.1, e = 0.2, N = 3000), all approaches yield quite comparable density profiles. 
Significant differences arise however in the death rates even at these large populations, 



c.f. Figure 10 



Our method applies to a whole class of models of branching random walks subject 
to a global constraint. A computational time-step of size e of such models consists of 
two sub-steps, which are very similar to the sub-steps of our simulation model. The first 
sub-step evolves the population density according to a branching random walk, which is 
regulated in the second sub-step to comply with a global constraint. Specifically, the first 
sub-step takes the concentration field c(x,t) to an intermediate state 

c(x, t + e) — c(x, t) = eCc{x, t) + y/2ec{x,t)ri(x, t) . (12) 

Here, the first term on the right-hand side describes the deterministic changes, which are 
due to net-growth and net-migration in our oasis model. In a time step e, the associated 
change in density is given by eC, where the Liouville operator C reads 

C = Dd 2 x + vd x + a6(x) - r . (13) 

Whereas the first two (particle conserving) terms describe biased diffusion with diffusivity 
D and average bias v, the last two reaction terms account for an overall negative growth 
rate — r except for the enhanced growth at the oasis. We note that, formally, we can 
reinterpret our model in an evolutionary context, if we consider x to be a coordinate 
measuring either a quantitative trait with continuous characteristics, such as height, 
skin, or body mass, or a co-ordinate measuring the genetic distance. In both cases, 
the biased diffusion corresponds to the random change due to mutations that often 
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tend to favor one direction. The delta-function growth spot would model an optimal 
phenotype which enjoys an increased fitness. The derealization threshold would have 
to be considered as Eigen's error threshold upon which selection is not able to maintain 
the optimal phenotype. More complicated quasi-species models can be easily devised, for 
instance, by considering a whole spectrum of mutations, recombination or more complex 
fitness landscape [18j H]. 



The stochastic second term on the left hand side in (12) accounts for the randomness 
in the reproduction process. The random function i](x,t) has vanishing mean and is 
uncorrelated across space and across the discrete time-steps, 



r)(x, t)r](x', f) = 5(x- x')5(t/ £ ),( t '/ £ ) • (14) 



The amplitude ^2ec(x, t) of the noise is chosen such that the produced variance represents 
the number fluctuations generated in the time span e by a population of size c(x, t), in 
which individuals branch and die at the same rate 1. The amplitude of the noise is thus 



consistent with the source of stochasticity in our simulations described in section [2TT 
In the second sub-step, 

c(x,t + e) = c(x,t + e)(l - A) , (15) 

the population density c(x, t + e) is rescaled by a factor 1 — A. Here, A has to be chosen in 
such a way that the density c(x, t + e) at the end of the time step satisfies the constraint 

dx u(x)c(x,t + e) . (16) 

The weighting function u(x) is at this point arbitrary. If one chooses u(x) = 1/N one 
obtains a fixed population size model, as we have simulated above (up to the lattice 
discreteness). For any other choice of u{x) the population size is not fixed, but becomes 
a fluctuating quantity. 



Notice that it is the constraint (16) that introduces a non-linearity in the problem, 
because A depends on c. The resulting noisy nonlinear reaction diffusion problems are 
intractable for general weighting functions u(x). This can be seen from the first moment 
equation for the mean c(x,t), for which one can derive the equation of motion 



d t c = (C - 2u) c - (c| (£t - 2u) u) c . (17) 



A detailed derivation of (17) can be found in [9]. Notice that the second term on the 



right hand side in general involves the second moment, c(x,t)c(y,t). The equation of 
motion for the second moment, on the other hand, involves the third moment, and so on. 
The resulting infinite hierarchy of moment equations is intractable unless one is willing 
to introduce an ad-hoc truncation scheme. 

However, it is possible to obtain a closed first moment equation for a special choice of 
the weighting function. Suppose, the solution u*(x) of the nonlinear differential equation 

(£ f - 2u,) u* = (18) 
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exists, and is taken as our choice of the weighting function. Then, the non-linear second 
term of the moment equation (17) vanishes identically. As a consequence, the dynamics 
of the mean is controlled by a solvable linear differential equation 

d t c = (C - 2u*) c . (19) 

This equation is similar to the deterministic equation except for the negative last term 
that accounts for number fluctuations. 

In summary, branching random walks occur whenever individuals proliferate, and 
are therefore a frequent element in both ecological and evolutionary models. A given 
ecological model often has an interpretation in an evolutionary context as well, if the 
dynamical variables are reinterpreted. It is usually difficult to analyze these models when 
populations are regulated in size, for instance by a logistic non-linearity. The algorithm 
described above provides a method to construct a solvable model of constrained branching 
random walks, in which the population size is fluctuating but finite. The algorithm to find 
and analyze this model consists of first solving equation (18) for the weighting function 
u*(x). Secondly, one has to determine the mean population density c(x,t) by solving 
the linear equation (19), which depends on the previously determined weighting function 
u*{x). It has been shown in [9] that this weighting function has a natural probabilistic 
interpretation: u*(x) is the fixation probability that the descendants of an individual 
sampled from position x will eventually take over the population. It has also been found 
in both Refs. [HI [8] that the tuned model is useful to determine those properties, which 
are only weakly dependent on the form of the constraint. In this way, it was possible to 
determine the speed of traveling waves in an primarily analytic way. In the following, we 
will use this approach of model tuning to derive the behavior of the death rate in our 
stochastic oasis model. 



3.3. Stochastic oasis model 



When the Liouville operator £ in (13) corresponding to our oasis problem is used, the 
steady state equations for u*(x) and for c st (x) take the form 



= Ddlu* — vd x u* + (a5(x) — r)u* — 2u 2 
= Dd 2 c + vd x c + (a8(x) — r)c — 2u*c . 



(20) 
(21) 



To simplify the notation in the following, we introduce rescaled variables. We use 
X = xa/D to measure length in units of D/a. We also introduce new variables 
for the population density C(X) = c(XD/a)a and the weighting function U*(X) = 
u*(XD / a)D / a 2 . The equations of motion 








dyU* - vd x U* + 



S(X)- 



d\c + vd x C + 



S(X)- 



p 



U, - 2U, 
C - 2UX 



2 

* ) 



(22) 
(23) 
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together with the constraint 

1 = I dX U*{X)C{X) (24) 



then form a closed set of equations that only depends on two control parameters, the 
rescaled velocity u = v/a and the rescaled death rate p = ADra~ 2 . [f] 

3.3.1. Stochastic oasis model without convection It turns out that the convection terms 



in (22) and (23) considerably complicate the analysis. Therefore, we first focus on the 
case without convection, v = 0, and treat the windy situation in approximation later on, 
in section 13.41 

The equation for U* with v — restricted to X > can be mapped onto the 
mechanical problem of a point particle in a potential. When we interpret X as time and 
U* as the location of the particle, its dynamics follows Newton's equation 

& X U* = -du«V(U*) , (25) 
where the potential V(Z7#) is given by 

V(U.) = -2U?/3-p(J?/8. (26) 
Using this analogy, we obtain the weighting function C/*(x) through the integral 

dX = - £ , (27) 

y/2(E-V(U.)) 

where E represents the sum of kinetic and potential energy of the moving particle. The 
value of E is fixed by the requirement that the "trajectory" (here proportional to the 
density profile) approaches zero infinitely slowly in the limit X — > oo. Thus, the total 
energy must equal the potential energy at [/* = 0. Hence, E = V(0) = 0. The integral 



in (27) from £/*(0) to U*(X) can then be performed analytically, which yields 

-2 



u,(x) = I 



sinh ^JpX/A — arcsch 



16C/*(0) 



(28) 



3p 

The value Z7* (0) has to be fixed by the boundary condition at X = 0. This boundary 



condition is obtained from the full equation for £/*, equation (22), by integrating from 
X = -e to X = e < 1, 

[Ul(e)-Ul(-e)} + U*(0) = 0, (29) 

which also holds for v > 0. In the absence of the convection term, we expect Z7* to be an 
even function. Therefore, we must demand the boundary condition 

U'M = -C/*(0)/2 . (30) 

| The rescaled death rate was chosen such that, in the deterministic limit, we have p = 1 in the convection 
less case. The rescaled velocity was chosen such that the deterministic derealization threshold is at v = 1. 
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(normalized) death rate p 



Figure 8: Mean population size N obtained for the analytical model with tuned 
constraint, in the absence of convection, v = 0. The population size diverges at 
p = 1, which corresponds to the deterministic limit. For p < p a w 0.44, the population 
profile c st acquires a pronounced algebraic decay up to a length £ ~ p~ x l 2 . This broad 
tail arises from large scale center-of-mass fluctuations. At p — p m — 1/16 the mean 
population size has a minimum. 



This boundary condition is satisfied by our solution (28) when we set 

O.(0) = ^(l-P) • 
Due to the symmetry in X, the full solution for U* is given by 

-2 



U m (X) 



3p 
16 



\fp 

sinh ( + artanh (>/p) 



(31) 



(32) 



Note that the form of the equations (22, 23) for £/*and C implies that U*(x) is proportional 
to a stationary density profile, 



C 



St 



11 



(33) 



where the proportionality constant (J dX U 2 (X)^j 1 is fixed by the constraint in (24). 

By integrating over the density profile, N = (D/a 2 ) J dX C st (X), we obtain a 
relation between the mean population size iV and the rescaled death rate p = ADra~ 2 , 

N=° 2 16 . (34) 
a 2 1 + 7p - 2p K J 

The predicted dependence of the mean population size on the death rate p is depicted in 
Figure [8] Notice that the (mean) population size N reaches a minimum at p = p m = 1/16. 
This minimum is a distinctive feature of the tuned model close to the derealization 
transition. In this regime, the population frequently detaches from the oasis and takes 
large excursions away from the oasis until it is eventually pushed back to the oasis. 
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During such excursions, the population explores the tail of the function u*. To obey 



the constraint (16) for a low value of it*, the population has to be rescaled to a large 
value, which leaves its footprint in an overall increased mean population size. Since this 
phenomenon cannot occur in a model with u =const., fixed population size models do 
not exhibit such a minimum in the population size. 

3.3.2. Deterministic and stochastic limits Our solutions for [/* and C s t are valid as 
long as 

ADr . . 

< p = —5- < 1 . 35 
cr 



The bounds in (35) can be readily interpreted. For p — >■ 1, we recover the exponential 



decay known from the deterministic limit (10) 



C st (X) ~exp(-|X|/2) . (36) 

The total population size scales as N ~ 32_D/(3a 2 5p) for Sp ee 1 — p <C 1. Since relative 
fluctuations in population size become small for large populations, we expect this leading 
order result to be independent of the form of the constraint. Indeed, the measured growth 
rates in the fixed N model for v = agree very well with our prediction, as can be seen 



from Figure 10 



In the case of very small death rates, p -C 1, we obtain an algebraic decay with 
distance, 

Cst(X) + 2 , (37) 

This algebraic decay indicates that the center-of-mass of the population fluctuates 
strongly around the center of the oasis. Figure [9] depicts the density profile for small 



but finite p, for which the decay is algebraic (as in (37)) up to a characteristic crossover 
scale (in original units) 

£~2W-. (38) 
V r 

For x ^> £, the decay is exponential. We note that to clearly see the algebraic tail, the 
rescaled growth rate has to be finely tuned to a very small value. Yet, an algebraic decay 
of correlations and a diverging correlation length, is reminiscent of a second order phase 
transition. This raises the question of whether the underlying physics is generic for a 
wide range of systems. 

Some more insight can be gained with the help of a simple scaling argument for 
the extinction dynamics close to the transition. For very small death rates, assume that 
the oasis sheds off subpopulations into the surrounding environment. Since this desert is 
only mildly deleterious, these trailing populations will survive up to a maximal time of 
roughly 1/r. Within this time, the subpopulations diffuse up to a length scale yjD/r. 



Incidentally, this is precisely the correlation length, defined in (38), where the density 
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profile crosses over from algebraic to a much faster exponential decay. The algebraic 
decay of the density profile within the correlation length can be explained by involving 
the statistics of branching random walks. Suppose, the populations are shed off from 
the oasis at a small, but constant rate k a ff. Each of these populations follows a critical 
branching process as long as the population is smaller then 1/r. (For larger populations, 
selection dominates and leads to a continual decrease in population size at rate r.) First, 
the probability for a small population to reach a size larger than N individuals through a 
critical branching process is proportional tol/N. The time to reach the population size N, 
conditional on survival, is roughly given by N generations. During this time, the position 
Xcom of the center-of-mass will have diffused by a length scale y/DN. Combining these 
scaling arguments, we determine the probability that a sub-population reaches further 
than the distance x as 

Pr(X com >x)~ Pr(T > x 2 1 D) ~ Pr(N > x 2 1 D) ~ ^ . (39) 

x 2 

Now, the mean population density at x will be given by the probability density to reach 
x, Pr(X com = x) ~ Dx~ 3 , times the typical density of a population that reaches x. 
Such a population has size N x ~ x 2 / D and a lateral spread of y/DN x (acquired during 
a coalescence time of N generations). Hence, the density of population that reaches x 
scales as x/D, and we obtain 

c(x) ~ Pr(X com = x)x/D ~ x~ 2 . (40) 

The pre-factor clearly depends on the rate at which populations are shed off from the 
oasis. 

The above scaling picture suggests, that the exponent —2 is generic for a critical 
branching random walk with a boundary condition of a finite (source) population on one 
end. In this context, critical means that birth and death rate are equal. If the landscape 
is only nearly critical and characterized by a weak selective imbalance r, then we expect 
to still see the power law up to the cutoff y^D/r. 



3.4- Stochastic oasis model with convection 

A finite wind speed significantly complicates the mathematical analysis, because it results 
in a non-hermitian Liouville operator, defined in (13). As long as the population remains 
localized, it is however possible to return to a hermitian problem by the non-linear variable 
transformation, 

U*(X) = e^X) , (41) 
C(X) = e-^$(X) . (42) 
In terms of the new functions and <&{X), the stationary equations of motion read 



= d 2 x V(X) + 
= d\<&{X) + 



5{X) 
5(X) 



p 


v 2 ' 


4 


4 


P 


v 2 ' 


4 


4 



2e^V 2 (X) , 
2e^(X)$(X) . 



(43) 
(44) 
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10 2 

position x 

Figure 9: Density profile c s t(x) close to the stochastic threshold, p < 1. The 

plot depicts the density profile in (33, 32) predicted from the solvable model with tuned 
constraint. We used the parameters a = 0.1, r = 10~ 7 ,D = 1, such that the rescaled 
death rate is much smaller than one, p = 4Dra~ 2 = 4 x 10 -5 . Notice the intermediate 
algebraic decay (with exponent —2) up to the correlation length £ ~ 2-i ' 



These equations cannot be solved in closed form, but a simple approximation is possible 
for small convection velocities where the exponential factors in ( 43|44 ) are close to 1. 



This approximation is justified if the exponent is smaller than 1 on the considered length 



scale. If we consider the localization length of the deterministic problem (c.f. (11)) as 
the relevant scale, this implies that the convection speed should be much less than the 
deterministic threshold, v <^ 1. 



After replacing the exponential factors in (43, 44) by 1, the equations of motion 
become identical to the convection free case with the replacement p t— > p+v 2 . Accordingly, 
the density profiles then obey 

C 8t {X;p,v) ^C st (X;p + v 2 ,0)e-^ , v < 1 . (45) 

The resulting predictions for the population densities and growth rates have been used 
in Figure [7] to compare with our fixed iV simulations. The agreement is not perfect, 
due to the combined error generated by the model difference and our small velocity 
approximation, but reproduces well the qualitative behavior of the simulations. 

3.4-1- Approximation for the death rate Finally, we use the tuned model to derive 
a rather simple prediction for the death rate p as a function of population size and 
convection speed. This approximation will help us better understand the derealization 
mechanism in the presence of number fluctuations. 
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Figure 10: Relation between death rate p and population size N for different 
convection velocities v. Symbols represent measured death rates in the fixed 
population size model for different population sizes and convection velocities (see legend). 
Dashed horizontal lines are the deterministic approximations ([9]), which strongly deviate 
from the simulations of finite populations. The solid red line is the predicted (mean) 



population size N for a given death rate p, c.f. (34), derived from the solvable model with 
tuned constrained in the convection-less case. Notice that these predictions also capture 
the quantitative trend of the simulation data for fixed population sizes (red symbols). 
The green and blue lines are approximations to the tuned model with convection, c.f. 
(48). Simulation parameters are a = 0.1, e = 0.2. 



Integrating the governing equation (21) for the mean density profile c(x) over the 
space coordinate x yields 



= — r J dx c(x) + ac(0) — 2 J dx u*(x)c(x) . (46) 
Since the first term is proportional to the total population size N, and the last term is 



fixed by the constraint (16), we obtain in general 



acO 2 

r = - = . (47) 

N N K J 

This result is identical to the mean field prediction except for the negative last term. 

Hence, the death rate is lower than the mean-field prediction by the term 2/N. As a 

consequence, smaller convection speeds than in the mean-field approximation are required 

for the death rates to cross 0, which signals the derealization of the population. 



For using equation (47) predict the death rates for a given parameter set, we still 



need to find the occupancy c(0) at the oasis. This task requires the non-trivial solution 



of the set of equations in (43). By employing the small- velocity approximation (45) to 
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express c(0) in (47), we obtain 




v 

Id ' 



(46 



This approximation (48) is used in Figures [5, 10 to explain fixed N simulation data, and 
ultimately also to predict the fuzzy phase boundary in Figure [2j as implicit condition 
when the death rates reaches zero. 



4. Discussion 



We have studied localization of a finite population in an ecological context. Within our 
model, individuals enjoy an effective growth rate benefit in one site of a linear chain of 
subpopulations. We asked the question under which conditions this preferred growth 
spot, which has been called "oasis", might be able to pin the population. We found that, 
under a model of globally regulated total population size, escape from the oasis inevitably 
occurs due to random number fluctuations. The escape time however sensitively depends 
on both the population size and the convection speed. The derealization dynamics may 
be summarized in the phase diagram of Figure [2] When the convection speed is small 
enough and the population large enough, escape is extremely unlikely and the population 
can be considered as effectively localized. Due to the exponential dependence of escape 
times on N and v , the cross-over from long to short derealization times is very sharp in 
the N-v plane. The measured density profiles converge to the exponential decay for large 
population sizes, as predicted from deterministic mean field theories [2j|3]. The agreement 
with the deterministic theory generally becomes poorer as escape events become more 
frequent, either through an increase in convection velocity or through a reduction in 
population size. Deviations are pronounced in particular in the expected death rate that 
has to be used to control the population. In the absence of convection, this death rate is 
found to continually decrease as the population size is lowered. This is in contrast to the 
population size independent death rate predicted by the mean-field theory. 

To analyze these findings, we have solved a variant of the model with a tuned 
population size constraint. In this model, a weighted sum over the occupancies of 
all demes is held strictly constant. The weighting function is chosen such that the 
first moment equation is linear. This model has the advantage of being analytically 
tractable. We found that the tuned model exhibited a sharp derealization transition for 
all population sizes. The phase transition line in the N-v plane could be determined in a 
simple analytical approximation, which is indicated by the dashed red line in the phase 
diagram Figure [2] Interestingly, we found that our tuned model exhibits a continuous 
derealization transition driven solely by number fluctuations, even in the absence of 
convection. Close to the transition, the density profile exhibits a power law tail with 
exponent —2 up to a correlation length that diverges at the transition. 
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It has to be remarked that the reason for the sharp phase transition in the tuned 
model is ultimately due to the inhomogeneity of the weighting function. Occupancy of the 
oasis is weighted more strongly than occupancy of lattice sites away from the oasis. This 
is in contrast to the fixed population size constraint, which has a spatially homogeneous 
weighting function: The occupancy of all demes are summed up, independent on their 
identity, and compared with the pre-fixed value. Therefore, the stochastic escape from 
the oasis is a characteristic effect of this degenerate weighting function. It remains to 
be seen, how one could explain analytically the exponential dependence of the typical 
escape time on population size and convection speed. Nevertheless, our simulations show 
that the narrow cross-over region from small to large escape times is located close to the 
sharp phase transition line of our tuned model in Figure [2} Furthermore, our tuned model 
predicts reasonably well the behavior of the death rate required to keep the population 
size finite. Compared with the mean-field approximation, the death rate is lowered by a 
value of roughly 2/N, in fair agreement with fixed N simulations, c.f. Figure [5] 

Our results also have a natural interpreted in an evolutionary context. Consider 
a population evolving in a fitness landscape characterized by an isolated fitness peak. 
Mutations cause random movements within this fitness landscape, typically biased in 
certain directions, and reproduction is often approximated by a branching process subject 
to the constraint of a fixed population size. The resulting quasi-species model of a finite 
population is a higher dimensional version of what we have considered in this paper, but 
otherwise analogous. We therefore expect our results could be useful for advancing the 
quasi-species theory, whose stochastic version has so far remained largely unexplored. 
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A. Different choices for the definition "population on the oasis" 



In our simulations, detailed in Section 2.1, the population always delocalizes from the 
growth spot if they are run long enough. Our analytically tractable model, on the other 
hand, describes populations that never fully detach from the oasis, provided the model 



parameter are within the regime in which the equations of motion, (22 21), have non 



trivial solutions. To compare both approaches, we needed to average observables only 
over those configurations in the simulation, which can reasonably be called "localized". 



In Figure Al, we present the death rate for several alternative choices for the definition 
of localized states. We termed states localized if the center-of-mass ("com") of the 
population was within a certain distance to the oasis and/or if the site of the oasis has 
an occupancy Co larger than a certain threshold value. The tested definitions gave very 
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similar results. For definiteness, we chose Cq > 1 for presenting averaged simulation data 
in the main text. Note that, with this definition, the death rate averaged over localized 
states cannot entirely decay to zero at the derealization transition. This together with 
the model differences are reasons for the discrepancy between fixed N simulations and 
our analysis of the tuned model, as seen in Figures [5 10 and Al 
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Figure Al: Different criteria for localization yield qualitatively similar results. 

In the main text, the population was defined to be localized at the oasis if at least one 
individual was present at the oasis, cq > 1. Different criteria are of course possible. Here, 
we compare the following conditions: "|com| < X": the center-of-mass (com) is within a 
distance X of the oasis. "(c_ 5 , c , c +5 ) > x": the occupancies in demes —5, (oasis) and 
5 exceeds x simultaneously. We also checked a combination of the last two conditions. 
Overall we find that our choice, en > 1, is in good agreement with the others, and the 
predictions from the tuned constraint model (black line). 



B. Relaxation spectrum in the case without convection 

In the main text, we have analyzed the steady state of the analytically tractable tuned 
model. Here, we discuss the relaxation towards the steady state for the case without 
convection. The relaxation dynamics is governed by the eigenmodes and eigenvalues of 
the linear differential equation for C [TJ [2J [20] ■ The following analysis shows that only 
a single bound state exists, together with a continuum of extended, delocalized states. 
These extended modes decay exponentially, with a rate faster than p/4 in (rescaled) time. 



The dynamics of the mean density is described by (c.f. (19) and (13)) 

d t C= (£-2<7*)C, (B.l) 
£ =d\-vd x + (5(A) - p/4 . (B.2) 

The mean density C can be decomposed into a basis consisting of (right) eigenfunctions of 
the linear operator (C — 2U*). Without convection this operator is hermitian, rendering 
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Figure Bl: Spectral properties of the oasis model, (a) Without convection, v = 0, 
the spectrum can be obtained analytically and consists of the single bound state at T = 
and a continuous line of eigenvalues, the open interval T G (—00, —p/4), for the extended 
delocalized states. In analogy with [2] we expect the convection to render the eigenvalues 
complex (in pairs). The real part of these eigenvalues is approaching the bound state 
with increasing convection v until the first extended state reaches the bound state at the 
the derealization velocity u c . (b) Bound state (r = 0) and one of the delocalized states 
(r = —p/2). Both states are normalized to fulfill Co(0) = C-p/2(0). 



all eigenvalues and eigenfunctions real valued. Here, we only treat the convection-less 
case, v = 0. The eigenfunctions {(r(X)} can then be chosen as an orthogonal basis, 
satisfying 

(£-2igCr = rCr, (B.3) 

J dx Cv{x)Cr{x) = if r ^ r . (b.4) 

In addition to (B.4) we can demand the normalization of (q, J dX(o(X) 2 = 1, which is 
used later in (B.8). 



By the spectral theorem we can write 
C(X,t) = J2^(T)e rt Cr(X) 



(B.5) 



characterizing the relaxation dynamics of perturbations in the mean density. In (B.5) 
the coefficients {/«n(r)} are the projections of the mean density C at initial time onto 
the eigenstates Cr, 



« (r)= / dX C(X,0)( T (X) . 



(B.6) 

Note that our constraint fixes the coefficient of the stationary eigenfunction with 
eigenvalue zero, 



dX U*{X)C{X,t) 



(B.7) 
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/3 5> (r)e r * I dX( (X)( r (X) 



0«o(O) . 



(B.8) 
(B.9) 



In going from (B.7) to (B.8) we used the fact that is proportional to Co with a 
proportionality constant (3. 

Inserting the operator in ( |B.2 ) into (B.3) leads to 



°x 



3p 
8 V 
P 



sinh 



X\ + artanh(y / p) 



Cr(X) 



(B.10) 



for X > and v — 0. 

The eigenfunction £r is a linear superposition, CrPO = Cr PO + ^Cr PO' of two 
functions and Cr > given by 



C< +) (X) oc 



1 fiF 

2+ — + (sinh (F(X)))- 2 



3p 



cosh 



4~ } 



(X)oc 




+ lcoth(r(X))sinh 



+ 1 Y(X) 



+ (sinh 



+ lcoth(F(X))cosh 




+ 1 Y(X) 



(B.n; 



sinh 



+ 1 Y(X) 



+ 1 Y(X) 



(B.12) 



Here, we have introduced the abbreviation Y(X) 
K is fixed by the two boundary conditions 
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X\ +artanh( v /p). The parameter 



X— >±oo 



(B.13) 

Cr(0) = -Cr(0)/2 . (B.14) 

Different cases for T have to be distinguished now. F > —p/4: and T < —p/4: render 
the term ^J^- + 1 in the eigenfunctions either real or purely imaginary, respectively. 
First, we treat the case V > —p/4. The boundary condition (B.13) is fulfilled only for 
K = —1. For the second boundary condition (B.14) we obtain T = as a necessary 
requirement, leading to CoPO oc [sinh(y(X))]~ 2 . This is the (already known) stationary 
state, Cgt{X) oc CqPO- F° r r < —p/4, the arguments of the hyperbolic functions 
are (purely) imaginary. Using the identity between hyperbolic functions with complex 
arguments and trigonometric functions, we obtain oscillating (but bound) functions Cr^) 
which automatically comply with the first boundary condition. Straightforward, but 
tedious, calculations yield the coefficient K, which depends only on V and p. 
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These calculations imply that the spectrum of the linear operator (C — is given 

by 

T G {0} U (-p/4, -oo) , (B.15) 
as illustrated in Figure |Bla T = is the single bound state, whereas we obtain a 



continuum of extended, oscillating functions, which correspond to the delocalized states. 

Analogous to previous results on the deterministic case [2], we expect that the real 
part of the eigenvalue of the extended states increases with increasing convection speed. 
Hence, perturbations relax more slowly than in the convection-less case. At the critical 
convection speed the relaxation time diverges. 
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